Determining molecule diffusion coefficients on surfaces from a locally fixed probe: 

On the analysis of signal fluctuations 
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Methods of determining surface diffusion coefficients of molecules from signal fluctuations of a 
locally fixed probe are revisited and refined. Particular emphasis is put on the influence of the 
molecule's extent. In addition to the formerly introduced autocorrelation method and residence time 
method, we develop a further method based on the distribution of intervals between successive peaks 
in the signal. The theoretical findings are applied to STM measurements of copper phthalocyanine 
(CuPc) molecules on the Ag(lOO) surface. We discuss advantages and disadvantages of each method 
and suggest a combination to obtain accurate results for diffusion coefficients. 
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I. INTRODUCTION 

To describe the growth kinetics of adsorbates on sur- 
faces, knowledge of diffusion coefficients of atoms and 
molecules on the surface is of vital importance. Different 
experimental techniques can be utilized to measure cor- 
responding diffusion coefficients. Two of the first meth- 
ods were Field Emission Microscopji^ and Field Ioniza- 
tion Microscopy)^ Another technique is the laser induced 
thermal desorption,4v£, where first an area of adparticles 
is cleared by desorption and then the successive repop- 
ulation measured by desorption. For scanning tunneling 
microscopy (STM) several modes of operation have been 
described: Taking scans of the surface with the adsorbed 
particles at different time instants to identify individual 
particle displacements,— recording current signals in time 
with a locally fixed STM tip^^ and coupling the tip to a 
single adsorbate and tracking its motions.^ Different pro- 
cedures to extract diffusion coefficients from these and 
other techniques have been reviewed by Barth in 2000.^ 

Analysis of island densities with predictions from rate 
equations for submonolayer growth kinetics is a further 
powerful method, which allows one also to deter- 
mine binding energies between adatoms or of adatoms 
to small clusters. An extension of this method to multi- 
component systems is given in Refs. [H andfisl 

In recent years the self-assembly of organic molecules 
on surfaces has attracted much attention in connection 
with the idea of developing electronic devices on a molec- 
ular levelji^ Molecules in organic surface growth often 
have high mobilities and sizes large compared to the lat- 
tice constant of the substrate. In this situation not all 
techniques are equally well suited. A convenient and 
powerful method is the recording of temporal current 
fluctuations with a fixed scanning tunneling microscope 
tipi^ or, in principle, of signal fluctuations from any 
locally flxed microscopic or other probe.-^ For example, 
for molecules on insulating substrates, the frequency or 
height fluctuations in an atomic force microscopy could 



be used. 

In order to extract diffusion coefficients from signal 
fluctuations of a probe, a theoretical basis for the anal- 
ysis of the measurement is necessary. One approach was 
developed for the auto-correlation function (ACF) of an 
STM signaliSiii and applied, on the basis of the corre- 
sponding power spectrum, to oxygen atom diffusion on 
Si(lll)fi^ and more recently to hydrogen atom diffusion 
on Si(lll).^^ Another approach was followed in Ref. [l9l 
by analyzing the distribution of peak widths in the signal, 
which correspond to residence times of single molecules 
under an STM tip. We will refer to it as the residence 
time distribution (RTD) method in the following. For de- 
termining absolute values of diffusion coefficients, a gen- 
eral problem is that the theoretical modeling involves a 
commonly not precisely known "detection function" with 
some finite range. However, when the extension of the 
diffusing molecules is larger than the detection range, this 
problem becomes essentially irrelevant and the size of the 
molecule enters as the relevant length scale. 

In this work we first revisit both the ACF and the 
RTD method. For the ACF method we study, as a 
new aspect, modifications implied by the finite size of 
molecules and for the RTD method we give a detailed 
account of the functional form of the RTD in different 
time regimes. In the RTD, time intervals are consid- 
ered where the molecule is under (or very close to) the 
tip. Corresponding time intervals are also included in 
the calculation of the ACF. It can therefore not be pre- 
cluded that the interaction of the molecule with the tip 
is influencing the results. A possible influence can be 
evaluated systematically by changing the bias voltage^i^ 
or the distance between tip and substrate. In order to 
ensure that the influence is small independent of adjust- 
ments, we develop a further method, which is based on 
the analysis of the distribution of peak-to-peak distances 
in the signal. This method is referred to as the interpeak 
time distribution (ITD) method in the following. 

All three methods are applied to diffusion measure- 
ments of copper phthalocyanine (CuPc) on the Ag(lOO) 
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surface previously reported in Ref. 'lij. We discuss ad- 
vantages and disadvantages of the three procedures men- 
tioned above and a possible combination of them, which 
enables us to improve the quality of results for the diffu- 
sion coefficient. 



II. AUTO-CORRELATION FUNCTION (ACF) 

We consider N molecules diffusing on a flat substrate 
area A and a probing tip of an STM placed at a fixed 
position above the surface. The lateral tip position 
is chosen as the origin of an a;-y-coordinate system on 
the surface. Each diffusing molecule a with its center 
at position ra{t) = {xa{t),ya{t)) gives a contribution 
ia{t) = G{ra{t)) to the total tunneling current I{t) = 
^a=i^ait), where G(r) is the "detection function". We 
regard the molecules as non-interacting, which should be 
a valid assumption at low concentrations c = N/A, or 
more precisely on time scales smaller than P /D, where 
I - is the mean intermolecular distance and c 

the number density of molecules. The autocorrelation 
function C{t) — {6I{0)SI{t)) of the current fluctuations 
SI{t) = I{t) — (/) can then be written as 



Accordingly, the detection function can be defined as "in- 
dicator function" of the molecule 



c{t) = Nmt{t)) 

= c J d\o J d2riG(ri)/(ri, t|ro)G(ro) , (1) 

where /(ri,t|ro) is the two-dimensional diffusion propa- 
gator 
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To provide a quantitative prediction for C{t), knowl- 
edge of the detection function G(r) is necessary. For the 
STM this requires a detailed treatment of the tunnel- 
ing problem, which is a long-standing problem since the 
first pioneering work by Tersoff and Hamann^S Based 
on a simple approach, Sumetskii and Kornyshe'\»i^ sug- 
gested a Gaussian detection function for small particles 
(adatoms) as a result of the tunneling barrier. For ex- 
tended molecules, the lateral structure of the electronic 
charge density needs to be taken into account. In fact, 
without referring to first-principle calculations for the 
specific system under consideration, it is difficult to ob- 
tain accurate expressions for the tunneling current. 

While such more detailed treatments are interesting 
in themselves, they are not really required for deter- 
mining the diffusion coefficient D. A more practical 
approach is to transform I(t) into a rectangular signal 
Ircc{t) — Q{I{t) ~ Ic), where 9(.) is the Heaviside jump 
function {&{x) — 1 for a; > and zero else), and Ic is a 
threshold current to exclude background noise. If we con- 
sider situations, where at most one molecule can be in the 
detection region under the tip, each single molecule con- 
tribution ia (t) to Jioc (t) also transforms to a rectangular 
signal i^it) e {0,1}, and /,ec(0 = Ec^Tit) & {0,1}- 



G(r) = 



if r e 
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where Ai is the set of center positions of the molecule 
which give rise to a tunneling current I > Ic- To get a de- 
scription that is independent of details of the molecule's 
shape, we assign a circle to the set M. The radius R 
of this circle can, for example, be determined by setting 
the covered area of the molecule equal to ttB?, by the 
gyration radius of the set Ai, or by some other reason- 
able requirement. The function G(r) from Eq. ([3]) is then 
approximated by 8(i? — |r|). 

The mutual exclusion under the tip implies that the 
molecules are no longer non-interacting as assumed when 
deriving Eq. ([T|). This leads to nonzero cross-correlations 
(*q'^(*)*/3''^(0)); a ^ p. While these cross-correlations 
could be treated approximately, we concentrate here on 
times smaller than the mean residence time ^ := 
E? /D of a molecule under the tip. In this regime we 
can set {iT {tVp" ^'^)) ^ for a ^ ^. Accordingly, the 
self-part N{i^^'^{t)i'^°^{Q)) appearing in Eq. ([T]) gives the 
correlation function Grcc(i) — (^roc(^)^rcc(O)) (i.e. with- 
out subtracting a term (/roc)^)- Evaluation of Eq. ([1]) 
then yields 
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where /q is the modified Bessel function of zeroth order. 

The knowledge of Groc(i) for times t < th, as pre- 
dicted by Eq. ^ for single-particle diffusion, is suffi- 
cient to determine D. For t > t^, collective proper- 
ties and the associated cross-correlations had to be taken 
into account. In a rough treatment one could factorize 
(«r(*)«™''(0)) - (^a'^WX^r^O))' "^^''^^ amounts to add 
{cttR^Y' to the expression in Eq. Q for t ^ tr. 



Application to CuPc on Ag(lOO) 

A relaxed structure of the CuPc molecule in vacuum 
after energy minimization, as obtained from a DFT cal- 
culation, is displayed in Fig. [U together with. A circle 
of radius R = 7.6 A is assigned to it. The relaxed struc- 
ture can be expected to be essentially unmodified when 
the CuPc molecule is adsorbed on the Ag(lOO) surface. 
STM studies combined with DFT calculations showed 
that transition metal phtalocyanines generally lie fiat on 
the Ag(lOO) surfacci^ This is in agreement with STM 
images of CuPc islands on Ag(lOO),^^ 

After deposition of CuPc molecules on Ag(lOO) up to a 
coverage of 10-15%, islands form and a quasi-stationary 
state is reached, where the rates of detachment from and 



FIG. 1. Relaxed structure of the CuPc molecule. Notice that 
the terminating bonds are connected to hydrogen atoms. The 
assigned detection area marked in gray is a circle with radius 
R = 7.6 A. It corresponds to the smallest circle that entails 
all nuclei. 



attachment to islands balance each other. In this quasi- 
stationary state the diffusing CuPc molecules have an 
effective small number density c. Then the tunneling cur- 
rent was measured for a stationary STM tip and constant 
bias voltage. The height of the tip was preset indirectly 
by the choice of setpoint current /setpoint and bias volt- 
age C/bias- Further details of the experiment are given in 
Ref. [11. 

The diffusing molecules cause fluctuations in the tun- 
neling current I{t) as shown in Fig. [2] (a) for a time in- 
terval of 2 ms at T = 192 K. A histogram of the tun- 
nelling current, determined for the whole time series of 
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FIG. 2. (a) Fluctuating STM tunneling current I{t) observed 
during CuPc diffusion on Ag(lOO) at T = 192 K (C/bias=1.7 V, 
/setpoint =20 pA). The dashed line marks the threshold value 
Ic, above which the distribution of the current values has 
no longer a Gaussian shape, see Fig. [3] (b) The rectangular 
signal /rcc(i) = Q{I{t) — Ic) associated with I{t). The widths 
of the rectangles give the residence times r and the distances 
between the rectangles give the interpeak times r'. 
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FIG. 3. Histogram of measured current values for T — 192 K 
after shifting the maximum to / = 0. The solid line marks 
a Gaussian fit with zero mean and standard deviation de- 
termined from the histogram for negative currents. The in- 
set shows an enlargement to demonstrate the determination 
of the threshold value Jc used for separating the diffusion- 
induced signal from the noise (see text). 

length 40 s, is depicted in Fig. [3l In this figure, the max- 
imum of the histogram was shifted to / = 0.^'^ Negative 
/ values are associated with noise and their distribution 
can be fitted to a half-sided Gaussian. Extending this 
Gaussian distribution to positive / values yields the solid 
line in Fig. [3l For small positive / values this curve fits 
very well the data, implying that these values can also be 
attributed to noise. The diffusing CuPc molecules cause 
deviations from the Gaussian distribution at larger / val- 
ues. To separate the diffusion-induced fluctuations from 
the noise, we define a threshold value Ic > 0, where the 
Gaussian histogram of the number of occurrences drops 
below one, cf. Fig. |31 This threshold is marked by the 
dashed line in Fig. [2][a) and the corresponding Ircc(t) is 
shown in Fig.[2I^b). Using this procedure we determined 
the ACF of /rcc(i) for seven different temperatures inves- 
tigated in the experiments.— 

Normalized ACFs Crec(i)/C'rcc(0) for two temperatures 
are shown in Fig. Hl^a) (symbols) and fits to them with 
Eq. (jl]) and under the constraint t < th are marked as 
solid lines. The resulting diffusion coefficients are dis- 
played in the Arrhenius plot in Fig.|4l[b) (symbols), where 
a least square fit of = Dq exp(— /Jq/ZcbT) to these data 
(solid line) yields an activation energy Ea — 30 meV and 
a pre-exponential factor Dq = 10"^ "' cm^/s. 



III. RESIDENCE TIME DISTRIBUTION (RTD) 

The residence times r are the time intervals between 
entrance and exit of a molecule into the detection region 
A4 of the probe, see Fig. [2] By sampling these r values, 
the RTD ^'(r) is obtained. 

For a theoretical description of the RTD we need to 
tackle the problem of the diffusion of a molecule center 
in a circle with radius R and absorbing boundaries, see 
Fig. E^. Since the molecule center cannot be placed at 
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FIG. 4. (a) Normalized ACFs Crec(t)/Crec(0) as a function of 
time for two temperatures (symbols) and best fits with Eq. Q 
for t < TR (solid lines), (b) Arrhenius plot of the resulting 
diffusion coefBcients (symbols) . A least square fit (solid lines) 
according to the Arrhenius law yield Ea ~ 30 meV and Do = 
10-«-^ cmVs. 



the absorbing boundary, its initial center position is con- 
sidered to be at a distance eR from the boundary. This 
means that the initial probability density of the center 
position is a delta-function concentrated on the ring with 
radius (1 — e)R. Physically the length cR may be associ- 
ated with an elementary step size of the molecule on the 
substrate. 

Due to the symmetric initial condition, the probabil- 
ity density for the molecule center position r at time t 
depends on r = |r| only and is given by (see Appendix) 



p{r,t) 



^ Jo(a„r/fl) Jo(an(l - e)) _„2j- 
^^^^ ^TT^e , (5) 

n— 1 



Ji(a„) 



where Jfc(.) are the Bessel functions of order k, and the 
an are the (positive) zeros of Jq, Jo{an) = with < 
ai < a2 < ■ ■ ■ 





FIG. 5. Illustration of the geometry used for the calculation 
of (a) the RTD and (b) the ITD. The dashed hues indicate 
the uniformly distributed initial probabilities on the rings dis- 
placed by eR (a) and (e'R) (b) from the absorbing boundary 
(solid lines). MC stands for the diffusing molecule center. 



With p{r,t) from Eq. (O, the calculation of the RTD 
amounts to calculating the first passage time distribu- 
tion with respect to the absorbing boundary. This can 
be carried out by applying standard techniques from ran- 
dom walk theory, as described, for example, in Ref. [13. 
The probability that the molecule center has not left 
the circle until time t is (/)(i) = 27r drrp{r,t). The 
probability that the center leaves the circle in the time 
interval [r, r -I- At] is (pir) — 4>{t + Ar), implying that 
the probability density ^{t) for the residence times t is 
iI){t) = —2TTdrjQ dr rp{r,T). This yields 



■0(7 



TR 



n=l 



Ji(an) 



(6) 



We note that the characteristic time (eR)^ / D sets a lower 
limit above which Eq. ([5]) becomes applicable, because 
one would have to refine the continuum treatment for 
smaller t. For the experiments analyzed here, this is not 
of relevance because such small t are not resolved. 

For T > (ei?)^/D the functional form of the solution 
changes with respect to the characteristic time tr the 
molecule needs to explore the circle area. For t ^ tr, 
we find 



iP{t) (X — ( — 

TR \TR 



-3/2 



(7) 



This power law can be understood by considering the 
(negative) time derivative of the probability that the 
molecule is next to the absorbing boundary at time t 
[i.e. the efflux rate, which equals iP{t)]. For r <^ tr this 
probability is given by the ratio of the explored bound- 
ary section (cx -^/r) to the explored area inside the circle 
(cx r), yielding ^{t) ^ -dr^/r/r r^^/^. 
For T ^ Tr, we find 



2 aiJo(ai(l-£)) dr/rn 



TR 



Ji(ai) 



(8) 
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This result can be interpreted by noting that for r ^ 
T/j, where the occupation probabihty of the molecule 
is spread over the circle, the efflux rate is essentially 
constant and given by the inverse of the time tj^ for a 
molecule to reach the boundary. A Poisson process with 
this constant rate yields ^'(t) ^ t^^ exp{— const, t /tr) 



Application to measurements 

For the diffusion of CuPc on Ag(lOO) the r values 
were extracted from the rectangular current signal as de- 
scribed in Fig. [2] Distributions ^'(r) of these r values 
are shown in Fig. IHJa) for two representative tempera- 
tures. By fitting the exponential tail with Eq. (|8]), we 
obtain the diffusion coefficient D. In a self-consistency 
check we have assured that the tail regime used for the 
fitting fulfills the requirement t ^ tji = /D. In addi- 
tion the value e can be calculated from the prefactor of 
the fit. We find e ~ 0.3, which means that eR is about 
the lattice constant a — 2.89 A of the Ag(lOO) substrate. 
Inserting D and e, Eq. ([6]) yields the full distribution 
iP{t), which is marked by the solid lines in Fig. |6l[a). 
The very good agreement of the theoretical prediction 
with the measured data demonstrates the reliability of 
the approach. 

The diffusion coefficient D is shown in Fig. [Hb) for 
the two temperatures in (a) together with the five further 
investigated temperatures in an Arrhenius plot. From the 
slope of the fitted line, we find an activation energy Ea = 
30 meV and a pre-exponential factor Dq = 10"^'^ cm^/s. 
These values agree with the ones of the ACF analysis. 

We note, however, that these values deviate from 
those reported in Ref. [HI {Ea — 81 meV and Dq = 
10~^'^ cm^/s). The difference is due to the fact that error 
bars were taken into account in the analysis performed in 
Ref. [19I. One point in the Arrhenius plot, which was de- 
termined from the data set measured for T = 222 K, had 
a particularly small error and was largely infiuencing the 
slope of the fit line because of its exposed position with 
respect to the other points. When excluding this par- 
ticular point from the fitting, values Ea — 38 meV and 
Dq — lO"^ *" cm^/s are obtained, in fair agreement with 
the present analysis. We have refrained from including 
error bars here because the fits of the RTD for different 
temperatures yielded comparable errors and small differ- 
ences between them seem to be insignificant with respect 
to other possible sources of errors, as, for example, mi- 
nor temperature fluctuations or the influence of spatial 
inhomogeneities in diffusion profiles, that are associated 
with the fact that islands act as sinks for the molecule 
diffusion. 



IV. INTER-PEAK TIME DISTRIBUTION (ITD) 

The inter-peak times r' are the time intervals between 
the end of a peak and the beginning of the next peak in 
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FIG. 6. (a) Residence time distribution for CuPc/Ag(100) 
at two temperatures (symbols). The solid line marks the re- 
sult predicted by Eq. ((5]) after determining the parameters D 
and e from the exponential tail region t ^ tr (cf. Eq. ^ 
and discussion in text), (b) Arrhenius plot of the extracted 
diffusion coefficients. A least square fit (solid line) with the 
Arrhenius law yields an activation energy Ea = 30 meV and 
a pre-exponential factor Do = 10"^ '' cm^/s. 



the rectangular signal, see Fig. [2l The statistics of them 
is, for small t', dominated by entrance and exit of the 
same molecule into the detection region A4. For calcu- 
lating the contribution of these return processes to the 
ITD we analyze the diffusion of a molecule center with 
initial distance e'i? from a circular absorbing boundary 
with radius R, see Fig. [5] For the probability density 
of the molecule center to be at position r at time t one 
obtains (see Appendix) 



, ^ /""^ da , , , 

P{r,t)^ / 2^ " [ar/R, a) 

Woia{l + e'),a) 



(9) 



J^{a)+Yi{a) 
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where Wo{x,y) = Jo(x)Fo(y) - JQ{y)yQ{x) and Yq is the 
zeroth order Bessel function of second kind. 

The ITD can be derived analogously to the treat- 
ment of the RTD by taking the time derivative of the 
integral of p(r, t) over the outer area with respect to 
the circle. In the present case it is more convenient to 
take the flow through the absorbing boundary, tp{T') = 
^ds • [—DWp{r,t)]r=R, which, when making use of the 
Wronskian [Jq{x)Yo{x) — Jo{x)Yq{x)] = 2/ttx, yields 



2 



da a 



Wo(a(l + eO,a) 
Ji{a)+Y,^{a) 



(10) 



For t' Tii, the asymptotic behavior for t' ^ is 



: exp 



At' 



-3/2 



(11) 

Accordingly, ip^r') rapidly rises for small r' and, after go- 
ing through a maximum -^max at t^^x = if'RY' I^D, ap- 
proaches a power law with exponent (-3/2). This power 
law has an analogous physical origin as the power law 
in the RTD, see the discussion in Sec. IIIII after Eq. ([T]). 
Here tr is the typical time where the molecule center in 
Fig. [S] realizes the finite extent of the detection area, or, 
in other words, where the molecules realizes its size. 
For t' Tr, Eq. PU)) can be approximated by 



V'(t') 



21n(l -l-e') 

T'Xn^iT'/TR) 



(12) 



The asymptotics ~ (rln^ t)'^ follows from the fact that 
for large r' the detection area becomes very small with 
respect to the area explored by the molecule. Accord- 
ingly, ^'(r') scales as the probability of first return to the 
origin of a two-dimensional random walk. 

The large r behavior predicted by Eq. ((T^ is, however, 
of limited use, because another molecule can enter the 
detection area before the molecule, which has left this 
area at last, returns to it. The memory to a molecule that 
leaves the detection area is lost on time scales of order 
P /D ~ 1/cD. On these time scales different molecules 
can be regarded as entering the detection area with a 
constant rate. This rate should scale with the inverse 
mean time D/l^ for a molecule outside the detection area 
to enter it. Hence in the limit of large r', an exponential 
distribution is expected, 

V'(t') ~ cD exp {-kttDct') , (13) 

where k is a constant of order unity. 



Application to measurements 

The inter-peak intervals were sampled from the cur- 
rent signal as explained in Fig. [5] and the ITDs for two 
temperatures are shown in Fig. [7Ka) . Note that for very 
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FIG. 7. (a) ITD for CuPc/Ag(100) for two temperatures. The 
inset shows the exponential decay at large times, where the 
straight solid lines (dashed lines in the main log-log plots) are 
fits with Eq. (|13p . yielding D (see text). Using these D values, 
the solid lines in the main plots mark the short time behavior 
after fitting with Eq. (|lip . (b) Arrhenius plot of the extracted 
diffusion coefficients for all seven investigated temperatures. 
A least square fit (solid line) with the Arrhenius law yields an 
activation energy Ea = 31 meV and a pre-exponential factor 
Do = 10-3-^ cmVs. 



small times the data indicate a plateau, which is not con- 
tained in Eq. (|lll) . A refinement of the continuum de- 
scription would be necessary to capture this behavior, 
analogous to the very small times in the RTD. For times 
T ^ Tr the power law predicted by Eq. (fTTj) is reflected 
by the straight line behavior in the double-logarithmic 
plot, before eventually the exponential decay according 
to Eq. ([T^ takes over. 

To obtain D from the ITD, we first concentrate on the 
long-time behavior. Setting k = 1, the rate cD of the 
exponential decay is provided by the slopes of the lines 
shown in the insets of Fig. [TJa). The concentration is 
determined from the probability cttR^ = '^^Ti/T for a 
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TABLE I. Strengths and weaknesses of the three methods for determining D. 



Characteristics 


ACF 


RTD 


ITD 


Signal processing 


® Convenient by FFT 


Peak widths affected by 7c 
value 


© Negligible influence of Ic 
on interpeak distances 


Tip influence 


Possible 


Possible 


© Less likely 


Assumptions 


Non-interacting particles 


None 


© None in short-time regime 
Non-interacting particles 
in long-time regime 



molecule to be in the detection area, where T = Ylii'^i + 
t[) is the total measurement time. For the seven data 
sets taken at different temperatures, we find c values in 
the range 0.6 — 5 x 10~^ A , corresponding to coverages 
e = 0.01 - 0.08% of the diffusing CuPc molecules. The 
resulting D values are shown in the Arrhenius plot in 
Fig.El^b) and yield an activation energy of Ea — 31 meV 
and a pre-exponential factor of Dq = 10^^'^ cm^/s. Note 
that the assumed value k = 1 affects only the prefactor 
Dq but not the activation energy Ea- 

With D obtained from the long-time behavior, we can 
fit the remaining part of the ITD with Eq. PT|) . Corre- 
sponding curves are shown as solid lines in the double- 
logarithmic plots of Fig. ^a) . They give a good agree- 
ment with the experimental data. The fits yield e' ~ 0.5, 
which is consistent with e ~ 0.3 obtained in the analysis 
of the RTD, see Sec. IIIII Accordingly, we find again that 
e'R is of the order of the lattice constant a, as one should 
expect. 

In principle, the part of the ITD dominated by the 
single molecule diffusion can also be used to determine D. 
For this we have to notice that the necessary refinements 
of the continuum theory to describe the behavior of the 
ITD left to the maximum 

3V6D ^ D 

'^'"^'^ " ^nil + e')eyHe'Rr ^ ' ^ ^ 

are not expected to yield larger values of the ITD. In fact, 
when considering jump dynamics of the molecules with 
a rate D/a^ for short times, the ITD should behave as 
~ {D /a^) exp{— const. Dr'/a^), i.e. the largest value of 
the ITD should be of order D/a^. Due to matching with 
the continuum treatment, we can identify V'max with the 
maxima seen in Fig[7Ja) . This then is a convenient way to 
determine D/{e'R)^, and knowing this value, to extract 
D by fitting the part right to the maximum predicted by 
Eq. (fTTI) . Application of this alternative method indeed 
yields values for D and e' in good agreement with those 
discussed above. 

V. COMPARISON OF THE METHODS 

For applications it is important to clarify how the three 
methods are best combined to obtain most accurate re- 
sults for D. To this end we need to evaluate the strengths 



and weaknesses of each method. Let us first note that 
all methods work with a single input parameter, which is 
the molecule radius R. All other variables arise from the 
analyses described in Sees. UlllIVI Because tb, — R^/D 
is actually determined, uncertainties with respect to R 
slightly affect the diffusion coefficient. 

Table U summarizes the advantages and disadvantages 
of the three methods. The ACF can be readily calculated 
by a fast Fourier transformation without caring about 
peak identification in the signal. A disadvantage is that 
only the short-time regime t ^ tji is governed by single 
particle diffusion, while an accurate theoretical descrip- 
tion of the crossover to the long-time regime, governed 
by collective particle diffusion, requires a careful consid- 
eration of the mutual exclusion of the molecules (and 
possibly other interaction effects) . The decay of the cor- 
relation function within this regime can be small, which 
then affects the accuracy of the D values resulting from 
the fitting. Another drawback is that the determination 
of the ACF includes time intervals where molecules are 
under the tip and possible interactions with the tip can 
thus have an influence on the diffusion properties. 

Both the RTD and ITD method require some more 
effort in preparation of the data, which is connected to 
the determination of the threshold current Ic and iden- 
tification of the peaks. Once Ic is set, both the peak 
widths and interpeak distances can be extracted simul- 
taneously. Note that any method of determining Ic is 
associated with some uncertainty. For the ITD this is 
no problem in practice, because the interpeak intervals 
are comparatively large. We have checked that, when 
taking the interpeak distances as time intervals between 
peak maxima, the results do not change significantly. For 
the RTD the uncertainty of Ic is a more severe problem. 
Because molecules diffuse slowly into the detection area, 
the peaks in the original tunneling current signal have 
rather flat flanks. As a result the peak widths change 
more sensitively with Ic than the interpeak intervals. 

An advantage of the RTD is that D can be determined 
solely by analyzing the exponential tail for large residence 
times. One should note, however, that it may be diffi- 
cult to obtain a good statistics in the tail regime, when 
the molecules are highly mobile or small. In this case 
the peaks are narrow and it could be difficult to resolve 
them accurately. In the RTD method the interaction 
with the tip can infiuence the residence times and in this 
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case one would not determine the free diffusion of the 
molecules on the substrate. By systematically changing 
the bias voltage, a possible influence can be reduced to a 
minimum. A strength of the RTD is that it is related 
to a single-particle problem. 

The ITD method has the advantage that tip-molecule 
interactions can be expected to have, if at all, a marginal 
influence on the interpeak times (irrespective of tip- 
substrate distances or bias voltages). For small interpeak 
times, the ITD is essentially also related to a single- 
particle problem. For large interpeak times, an approx- 
imate value of D can be obtained based on an estimate 
for processes, where one molecule in the detection area 
under the tip is followed by another molecule. The values 
obtained in this way turn out to be close to those result- 
ing from the other methods. Compared to the residence 
times, the interpeak times are quite large and are thus 
less prone to the experimental time resolution and the 
threshold value Ic- 
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FIG. 8. Comparison of the diffusion coefficients obtained 
from the three methods. Squares, dashed line: ACF method 



{Ea = 30 meV, Do 



lo- 



om /s); Circles, dotted line: 



RTD method {Ea = 30 meV, Do = 10"^" cm^/s); Dia- 
monds, dashed-dotted line: ITD method {Ea = 31 meV, 



Do 



lo- 



om /s). 



VI. CONCLUSIONS 

In order to exploit the strength of each method one 
can combine the different methods, if the corresponding 
data in the relevant time regime have sufficient statis- 
tics. For a first classification of the measurement it is 
helpful to choose the ACF method, because it yields fast 
results without analyzing the time series in detail. In 
a subsequent step the peak widths and interpeak inter- 
vals should be extracted with the procedure described 
in Sec. |lll We then recommend to use the RTD method 
for determining D whenever the exponential tail is suffi- 
ciently pronounced and the tip influence on the diffusion 
can be neglected. Otherwise the ITD method should be 
preferred. In any case, one should always apply both 
the RTD and the ITD method to perform consistency 
checks and to obtain most reliable results. Our analysis 
for CuPc on Ag(lOO) shows that, despite diffusion coef- 
ficients extracted from the different methods may differ 
for one or another sample, the activation energies deter- 
mined from the Arrhcnius plots should have comparable 
values, see Fig. |S1 

An interesting question is whether the methods dis- 
cussed here can be taken over to other fields. Obviously, 
the ACF method is well known in fluorescence correla- 
tion spectroscopy (FCS)r^ where intensity fluctuations 
reflect concentration fluctuations, typically in some fi- 
nite detection volume. Similarly, analysis of density fluc- 
tuations has been applied as one variant to determine 
diffusion coefficients in Field Emission Microscopy.^ In 
the situation discussed in Sec. |TT] only one particle is in 
the detection area and accordingly information on the 
tagged particle diffusion (tracer diffusion) is obtained. In 
the hypothetical limit of very small particle concentra- 
tions in FCS, where the mean interparticle distance be- 
comes larger than the linear size of the detection volume, 
one should essentially recover the behavior discussed in 



Sec. |TT] (with R then playing the role of the detection 
length). In common applications of this method one is, 
however, not able to extract the tagged particle infor- 
mation and this makes a difference to the ACF method 
discussed in Sec. [Til This also prevents the use of the 
RTD and ITD in the analysis of typical fluorescence sig- 
nals. These methods can yet become useful in light of 
the ongoing development of sophisticated techniques to 
measure single-molecule properties. 

In summary, we provided formulas for different means 
of analyzing fluctuations in the signal of a locally fixed 
probe on a surface, where the fluctuations are caused 
by diffusing particles. By exploiting the fact that the 
molecule extension enters as a scale bar into the equa- 
tions, a quantitative determination of the diffusion con- 
stant becomes possible. So far, we have considered cir- 
cular shapes for planar molecules. The situation will ev- 
idently be more involved when rectangular or even more 
complicated shapes are considered. In addition rotational 
diffusion and anisotropics induced by the substrate then 
can modify the results. As mentioned above, the inter- 
action of the STM tip with the diffusing molecule was 
neglected so far. The interaction is influenced by the ex- 
perimental setup (tip shape, tip-substrate distance etc.) 
and it is specific for the tip and molecule material. One 
objective could be to minimize the tip influence exper- 
imentally. On the other hand, it could be interesting 
to study how the interaction potential affects the diffu- 
sion of the molecule. This potential has so far not been 
included in the analysis, but it should be possible and 
bears an interesting route to obtain information on the 
tip-molecule interaction. 
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Appendix: Diffusion propagators for RTD and ITD 

The results given in Eqs. ([5]) and ^ for p{r,t) have 
been derived earUer in the hterature, in particular Eq. © 
by using the Heaviside methodj2£ We give a straightfor- 
ward derivation here, using separation of variables and 
eigenfunction expansions. 

Let us first consider the equation for diffusion of a par- 
ticle in a circular stripe a < r < c {a here not equal to 
the lattice constant used in the main text) with absorbing 
boundaries and initial distribution 



Pair) 



1 

27rr 



S{r ~b) , a <b < c. 



(A.l) 



Due to the radial symmetry p{r,t) — p(r, t)/27r, where 
p{r, t) satisfies the radial diffusion equation 



dp{r,t) 
dt 



with p{a,t) — p{c,t) = and p(r, 0) = S{r — b)/r. 

Making the product ansatz p{r,t) — p{r,0)g{t), 
Eq. (|A.2p separates in the variables r and t. One obtains 
g{t) = where the allowed values for are the 

eigenvalues of the radial Laplace operator in the circular 
stripe, 

(A.3) 

Because the Laplacian is negative definite, > 0. 
The eigenfunctions are given by linear combinations 
of the zeroth order Bessel functions Jo(-) and Yo{.) of 
first and second kind, 4'„(r) = yl„ Jo(A„r) -I- BnYo{Xnr) 
(where restriction to A„ > gives linear indepen- 
dent eigenfunctions). The Dirichlet boundary condi- 
tions yield A„Jo(A„a) = -BnYo{Xna) [or A„Jo(A„c) = 
-i?„Fo(A„c)], and 



Jo(A„a)Fo(A„c) - Jo(A„c)ro(A„a) 







(A.4) 



as determining equation for the A„, n = 1,2,... (0 < 
Ai < A2 < . . .). Introducing the cross product 



W{x,y) = Mx)Yoiy) - Jo{y)Yo{x) , 
the solution becomes 



p{r, t) C„Wo(A„r, A„a)e" 



-A? Dt 



(A.5) 



(A.' 



where Wo(A„a, A„c) = 0. Utilizing the orthogonality of 
the eigenfunctions. 



dr r Wo{X„ir, X„ia)Wo{Xnr, A„a) 

= 6,nn I dr r tVg^(Anr, A„a) , (A. 7) 



the expansion coefficients C„ follow from the initial con- 
dition Eq. (|XT|) : 



dr r Wo{X„r, An, a)p(r, 0) 
J^dr r W§{Xnr,Xna) 
Wo(A„5, A„a) 



J^dr r W§{Xnr, A„a) 
The result for p{r, t) thus is 



(A., 



p{v,t) ^ y ^o(A»^A„a)W-o(A„r,A„a) 

h '^^SldrrW^iXr.r^Xr.a) ^ ' 

The diffusion propagator in Eq. ^ relevant for the ITD 
corresponds to the limit c — > 00, where the spectrum of 
eigenvalues determined by Eq. (IA.4I) becomes continuous. 
Analogous to the change of a Fourier series to a Fourier 
integral, we are led to consider the Weber transform^^ of 
a function f = f(r) 



F(A) = ^ /" drrWo{Xr,Xa)f{r) 

^ J a 

with back-transformation 

Wo(Ar, Xa) 



(A.IO) 



fir) = a' 



dXX- 



JO JiiXa)+YoH\a) 
Accordingly, Eq. (|A.6[) becomes 



^^(A). (A.ll) 



p{r,t)^ I dAC(A)I¥o(Ar,Aa)e 
"'0 



-VDt 



(A.12) 



where 
C(A) 

This yields 

Pir.t) = 



A 



JliXa)+YgiXa) A 

Aiyo(A6, Aa) 
JliXa)+Y^iXa) 



drrWa (Ar, Aa)p(r, 0) 

(A.13) 



dA Wo(Ar,Aa)Wo(Ab,Aa) _;,2^, 

2^ 4iXa)+YiiXa) 

(A.14) 

Equation © follows by setting a ^ R, b ~ il + e)R, and 
A = a/R. 

For the diffusion propagator in Eq. ([5]) relevant for 
the RTD only one boundary condition p(c, t) — has 
to be taken into account. In this case the Bessel func- 
tions of second kind cease to apply, because their log- 
arithmic singularity at the origin eliminates them from 
the space of functions, where the radial Laplace oper- 
ator is Hermitean. Notice that the logarithmic singu- 
larity would be no problem with respect of the integra- 
bility of p(r, t) . The eigenfunctions thus are given by 
^nif) = AnJoi^nf), whcrc the A„ are determined by 
Joi^nc) = 0. Equation (|A.6p becomes 



p(r, t) ^y A„ Jo(A„r) e" 



KDt 



(A.15) 
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and the An are again determined by the initial condition, This yields 
corresponding to an expansion of p{r, 0) into a Fourier- 

Bessel series, , ^ J_ Jo(Anr) Jo(A„6) ^-a^d* ^j-. 

Trr2 Ji(A„c) 



n=l 



Equation ([5]) follows by setting c — R, b — {1 — e)R and 



/ ;dr r Jo(A„r)p(r,0) ^ 2Jo(A„6) 
/;drrj2(A„r) c2j2(A„c) 
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